#pragma  warning (disable:4244) // loss of precision when converting double to float
#pragma  warning (disable:4305) // truncation from 'double' to 'const float'

#define _USE_MATH_DEFINES
#define PRECISION float

template<typename T>
inline T sq(T x) {return x*x;}

/*
OpenMesh::Vec3f x,y,n,crossproductXY;
...
l = (x-y).length();
n = x.normalize();
scalarProductXY = (x | y);
crossProductXY = x % y;
*/

#include <iostream>
#include <levmar.h>
#include <OpenMesh/Core/IO/MeshIO.hh>
#include <OpenMesh/Core/Mesh/PolyMesh_ArrayKernelT.hh>
#include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
#include <OpenMesh/Tools/Subdivider/Uniform/LoopT.hh>
#include <OpenMesh/Core/Mesh/Status.hh>
#include "polygon.h"
#include "triangle.h"
#include "vector.h"

const float PITOGRAD = 180 /(float)M_PI;


// ----------------------------------------------------------------------------

struct MyTraits : public OpenMesh::DefaultTraits
{
  VertexAttributes(OpenMesh::Attributes::Status);
  FaceAttributes(OpenMesh::Attributes::Status);
  EdgeAttributes(OpenMesh::Attributes::Status);
};


#include <Imath/ImathMatrix.h>
#define IM IMATH_INTERNAL_NAMESPACE
typedef IM::Matrix33<PRECISION> Mat3;

// ----------------------------------------------------------------------------

using namespace std;
typedef OpenMesh::TriMesh_ArrayKernelT<MyTraits> TriMesh;
typedef OpenMesh::PolyMesh_ArrayKernelT<>  PolyMesh;

Mat3 pointsToMat(OpenMesh::Vec3f& A, OpenMesh::Vec3f& B, OpenMesh::Vec3f& C)
{
	Mat3 m;
	m[0][0] = A[0]; 	m[0][1] = A[1]; 	m[0][2] = A[2];
	m[1][0] = B[0]; 	m[1][1] = B[1]; 	m[1][2] = B[2];
	m[2][0] = C[0]; 	m[2][1] = C[1]; 	m[2][2] = C[2];
	return m;
}


bool readmesh(TriMesh&mesh, char* filename)
{
	return OpenMesh::IO::read_mesh(mesh, filename) ;
}

void subdivide(TriMesh& mesh, int k)
{

	// Initialize subdivider
	OpenMesh::Subdivider::Uniform::LoopT<TriMesh> loop;
	loop.attach(mesh);
	loop( k , false);
	loop.detach();
} // void subdivide()

double perimeter(OpenMesh::Vec3f& pointA, OpenMesh::Vec3f& pointB, OpenMesh::Vec3f& pointC)
{
	double d;
	d = (pointB-pointA).length();
	d += (pointC-pointB).length();
	d += (pointA-pointC).length();
	return d;
}

double area(OpenMesh::Vec3f& pointA, OpenMesh::Vec3f& pointB, OpenMesh::Vec3f& pointC)
{
	double d;
	d = 0.5 * cross(pointB-pointA , pointC-pointA).norm();
	return d;
}

double volume(OpenMesh::Vec3f& pointA, OpenMesh::Vec3f& pointB, OpenMesh::Vec3f& pointC)
{
	// http://www.ditutor.com/vectors/volume_tetrahedron.html
	Mat3 m = pointsToMat(pointA , pointB, pointC);
	// use minus sign or change the order of the points in above function call 
	return -m.determinant()/6.0;
}


void analyzeTriangle(OpenMesh::FaceHandle & _fh, TriMesh& mesh)
{
	OpenMesh::Vec3f pointA , pointB , pointC;
	TriMesh::ConstFaceVertexIter cfvIt;
	cfvIt = mesh.cfv_iter(_fh);
	pointA = mesh.point(cfvIt.handle());
	pointB = mesh.point((++cfvIt).handle());
	pointC = mesh.point((++cfvIt).handle());
//	cout<<perimeter(pointA, pointB, pointC) << "   "<< area(pointA, pointB, pointC)<<endl;
	cout<<volume(pointA, pointB, pointC)<<endl; // volume of the tetrahedron with vertices A, B, C and the origin (0,0,0)
	
	//cout<<pointA<<endl;
}



float angle(OpenMesh::Vec3f& pointA, OpenMesh::Vec3f& pointB, OpenMesh::Vec3f& pointC)
{
	float sinTheta = cross(pointA-pointB, pointC-pointB).norm() / ((pointA-pointB).norm()*(pointC-pointB).norm());
	float cosTheta = dot(pointA-pointB, pointC-pointB);
	float theta = asin(sinTheta)*PITOGRAD; 
	if (cosTheta<0)
	{
		theta = 360 - theta;
	}
	return theta;
}

double min3(double a, double b, double c)
{
	if (a<b) return min(a,c); else return min(b,c);
}



TriMesh::FaceHandle addFaceSimple(TriMesh& mesh, TriMesh::VertexHandle a, TriMesh::VertexHandle b, TriMesh::VertexHandle c)
// adds a new face to the mesh, if the verteces are already part of the mesh
{
	TriMesh::FaceHandle fh = mesh.add_face(a, b,c );
	if (fh.idx()<0) // invalid face
		fh = mesh.add_face(b, a, c); // reorder the points
	return fh;
} // addFaceSimple


TriMesh::FaceHandle addFace(TriMesh& mesh, TriMesh::VertexHandle a, TriMesh::VertexHandle b, TriMesh::VertexHandle c)
// adds a new face to the mesh, if the verteces are already part of the mesh
{
	static std::vector<TriMesh::VertexHandle>  face_vhandles(3);
// compute quality
	OpenMesh::Vec3f pointA = mesh.point(a);
	OpenMesh::Vec3f pointB = mesh.point(b);
	OpenMesh::Vec3f pointC = mesh.point(c);

	double angleA = angle(pointC, pointA, pointB);
	double angleB = angle(pointA, pointB, pointC);
	double angleC = angle(pointB, pointC, pointA);
	float ab = (pointA-pointB).norm();
	float bc = (pointB-pointC).norm();
	float ca = (pointC-pointA).norm();

	OpenMesh::Vec3f av = pointA-pointC;
	OpenMesh::Vec3f bv = pointB-pointC;
	// http://en.wikipedia.org/wiki/Circumcenter
	float circumdiameter = av.norm()*bv.norm()*(av-bv).norm() / cross(av,bv).norm();
//	double quality = min3(angleA, angleB, angleC);
	double quality = min3(ab, bc, ca) / circumdiameter; // bad if quality < Threshold
	cout<< "Adding" << a <<"  " << b << "  "<< c<<"  "<<quality<<endl;

//	if (quality<threshold) addFaceToListofFutureImprovements;
	return addFaceSimple(mesh, a,b,c);
} // addFace


// deletes the triangle identified by <f> and adds a new tetrahedron to the mesh formed by the vertices of the deleted tringle and the new point <P>
void extendMesh(TriMesh& mesh, TriMesh::FaceHandle& fh, OpenMesh::Vec3f& P, bool runGarbageCollector=true)
{

	OpenMesh::VertexHandle newP = mesh.add_vertex(P);
	TriMesh::ConstFaceVertexIter cfvIt;
	cfvIt = mesh.cfv_iter(fh);
	OpenMesh::VertexHandle pointA = cfvIt.handle();
	OpenMesh::VertexHandle pointB = (++cfvIt).handle();
	OpenMesh::VertexHandle pointC = (++cfvIt).handle();
	mesh.delete_face(fh, false);
	addFace(mesh, pointA, pointB, newP);
	addFace(mesh, pointB, pointC, newP);
	addFace(mesh, pointC, pointA, newP);
	if (runGarbageCollector)
	{
		mesh.garbage_collection();
	}
}

TriMesh createTetra(double edgelen)
{
	TriMesh mesh;
	TriMesh::VertexHandle vhandle[4];
	vhandle[0] = mesh.add_vertex(PolyMesh::Point(edgelen/2, 0,  -edgelen*sqrt(2.0)/4.0));
	vhandle[1] = mesh.add_vertex(PolyMesh::Point(-edgelen/2, 0,  -edgelen*sqrt(2.0)/4.0));
	vhandle[2] = mesh.add_vertex(PolyMesh::Point( 0,  -edgelen/2,  edgelen*sqrt(2.0)/4.0));
	vhandle[3] = mesh.add_vertex(PolyMesh::Point(0,  edgelen/2,  edgelen*sqrt(2.0)/4.0));

	//std::vector<TriMesh::VertexHandle>  face_vhandles;

	addFace(mesh, vhandle[1], vhandle[2], vhandle[0]);
	addFace(mesh, vhandle[0], vhandle[2], vhandle[3]);
	addFace(mesh, vhandle[3], vhandle[2], vhandle[1]);
	addFace(mesh, vhandle[1], vhandle[0], vhandle[3]);

	return mesh;
}// TriMesh createTetra()

#include "tutorials.cpp"

double meshVolume(TriMesh& mesh)
{
	mesh.request_face_normals();
	mesh.update_normals();
	double vol = 0;
	double surf = 0;

	//iterate through all faces;
	TriMesh::FaceIter f_it, f_end(mesh.faces_end());
	OpenMesh::Vec3f pointA , pointB , pointC;
	TriMesh::ConstFaceVertexIter cfvIt;

	for (f_it = mesh.faces_begin(); f_it!=f_end; ++f_it)
	{
		cfvIt = mesh.cfv_iter(f_it);
		pointA = mesh.point(cfvIt.handle());
		pointB = mesh.point((++cfvIt).handle());
		pointC = mesh.point((++cfvIt).handle());
		surf += area(pointA, pointB, pointC);
		vol += volume(pointA, pointB, pointC);
	}
	return vol;

} // double meshVolume(TriMesh& mesh)


OpenMesh::FaceHandle faceClosestToPoint(TriMesh& mesh, OpenMesh::Vec3f& P)
{
// return a handle to the face that is closes to a give point.
// distance is measured as sum of Euclidean distances from the 3 vertices to the point.
	vector<float> dist(mesh.n_vertices());
	unsigned int i = 0;

	TriMesh::VertexIter v_it, v_end(mesh.vertices_end());
	for (v_it = mesh.vertices_begin(); v_it!=v_end; ++v_it)
	{
		dist[v_it->idx()]=(P-mesh.point(v_it)).length();
	}

	TriMesh::FaceIter f_it, f_end(mesh.faces_end());
	TriMesh::ConstFaceVertexIter cfvIt;
	double d, mind=10000000000;
	OpenMesh::FaceHandle mindHandle;
	for (f_it = mesh.faces_begin(); f_it!=f_end; ++f_it)
	{
		cfvIt = mesh.cfv_iter(f_it);
		d = dist[ cfvIt.handle().idx()];
		d += dist[ (++cfvIt).handle().idx()];
		d += dist[ (++cfvIt).handle().idx()];
		if (d<mind)
		{
			mind = d;
			mindHandle = f_it.handle();
		}
	} //for (f_it = mesh.faces_begin(); f_it!=f_end; ++f_it)
	return mindHandle ;
} // OpenMesh::FaceHandle faceClosestToPoint(TriMesh& mesh, OpenMesh::Vec3f& P)


TriMesh createSphere(float radius, int subdivisionLevel)
// creates a sphere centred in (0,0,0) with radius "radius" 
// See http://en.wikipedia.org/wiki/Tetrahedron for radius of circumsphere formula
{
	float tetraEdge = 4 * radius / sqrt(6.0);
	TriMesh mesh = createTetra(tetraEdge); // OpenMesh::IO::write_mesh(mesh, "tetra.off"); readmesh(mesh, "tetra.off");
	subdivide(mesh, subdivisionLevel);
	OpenMesh::Vec3f O(0,0,0);
// iterate through all vertices and push them to the sphere
	TriMesh::VertexIter v_it, v_end(mesh.vertices_end());
	OpenMesh::Vec3f P;
	for (v_it = mesh.vertices_begin(); v_it!=v_end; ++v_it)
	{
        P = mesh.point(v_it);
		P *= radius / P.norm();
		//cout<<P<<" --> "<<P.norm()<<endl;
		mesh.set_point(v_it, P);
	}
	return mesh;
}


typedef LightVector<TriMesh::HalfedgeHandle> hehV;

void findHoles(TriMesh& mesh, LightVector<hehV>& holes)
{
	// TODO delete isolated vertices to be sure


	// iterate through all edges
	unsigned int n = mesh.n_halfedges();
	unsigned int i;
	int* edgeinfo = new int[n];

	for (i=0; i<n; ++i) 
	{
		edgeinfo[i]=-1;
	}

	TriMesh::HalfedgeIter he_it, he_end(mesh.halfedges_end());
	TriMesh::Halfedge H;
	TriMesh::HalfedgeHandle heh, heh_init;

	hehV hole;
	int nbHoles(0);

	for (he_it = mesh.halfedges_begin(); he_it!=he_end; ++he_it)
	{
		heh_init = he_it.handle();	
		H = mesh.halfedge(heh_init);
		i = heh_init.idx();
		if (! mesh.is_valid_handle(H.face_handle_)) // boundary
		{
			if (edgeinfo[i]==-1)
			{
				//cout<<"Hole ["<<nbHoles<<"] = "<<heh_init;
				hole.push_back(heh_init);
				edgeinfo[i] = nbHoles;
				// We can do this as often as we want:
				heh = mesh.next_halfedge_handle(heh_init);
				while(heh != heh_init) // go through the boundary edges
				{
					//cout<<" "<<heh;
					hole.push_back(heh);
					edgeinfo[heh.idx()]=nbHoles;
					heh = mesh.next_halfedge_handle(heh);
				}
				++nbHoles;
				//cout<<endl;
				holes.push_back(hole);
				hole.clear(); // get ready for next hole
			}
		} else
		{
			edgeinfo[i]=-1; // this is not a boundary edge
		} // if (! mesh.is_valid_handle(H.face_handle_)) // boundary
	} // for (he_it = mesh.halfedges_begin(); he_it!=he_end; ++he_it)
	delete[] edgeinfo;
} // void findHoles(TriMesh& mesh)


void displayHoles(TriMesh& mesh, LightVector<hehV>& holes)
{
	unsigned int i, j;
	for(i=0; i<holes.size(); ++i)
	{
		cout<<"Hole "<< i<<endl;
		for(j=0; j<holes[i].size(); ++j)
		{
			cout<<holes[i][j]<<" "<<mesh.point(mesh.halfedge(holes[i][j]).vertex_handle_)<<endl;
		}
		cout<<endl;
	}
} // void displayHoles(vector<hehV>& holes)


void printOff(TriMesh& mesh, hehV& hole)
{
	//std::ostream& io = cout;
	ofstream io;
	io.open("outHole.off");
	unsigned int i;

	io<<"OFF"<<std::endl;
	io<< hole.size() << "  " << hole.size() << "   0" << std::endl;

	OpenMesh::Vec3f P;
	for (i=0; i<hole.size(); ++i)
	{
		P = mesh.point(mesh.halfedge(hole[i]).vertex_handle_);
		io<<P[0]<<"  "<<P[1]<<"  "<<P[2]<<endl;
	}
	// writing degenerate triangles as edges
	for (i=0; i<hole.size()-1; ++i)
	{
		io<<"3 "<<i<<" "<<i+1<<" "<< i<< endl;
	}
	io<<"3 "<<0<<" "<<hole.size()-1<<" "<< 0<< endl;

	io.close();
}




double localDensity(TriMesh& mesh, TriMesh::VertexHandle& vh)
{
	// compute localDensity as average distance between the point and its 1-ring neighbourhood

	int valence = 0;
	double dist = 0;
	TriMesh::Point P, O;
	
	O = mesh.point( vh );
	
	TriMesh::VertexVertexIter    vv_it;
	for (vv_it=mesh.vv_iter( vh ); vv_it; ++vv_it)
	{
		P = mesh.point( vv_it );
		dist += (P-O).norm();
		++valence;
	}
	return dist/valence;
} // double localDensity(TriMesh::VertexHandle& vh)



double angleBetweenPoints(TriMesh& mesh, TriMesh::VertexHandle& Ah, TriMesh::VertexHandle& Oh, TriMesh::VertexHandle& Bh)
{
	// returns the angle formed by the points AOC
	TriMesh::Point A, O, B;
	A = mesh.point( Ah );
	O = mesh.point( Oh );
	B = mesh.point( Bh );
	double a = (A-O) | (B-O);
	return a;
}

bool closeHole(TriMesh& mesh, hehV& hole)
{
//	TriMesh::VertexHandle vhandle[3];
	unsigned int n = hole.size();
	unsigned int i;

	if (n<3) throw "Big ugly error.";

	// only one triangle missing
	if (n==3)
	{
		addFace(mesh, mesh.halfedge(hole[0]).vertex_handle_, mesh.halfedge(hole[1]).vertex_handle_, mesh.halfedge(hole[2]).vertex_handle_);
		return true;
	}
// just zip it up - trivial ugly solution
//goto alternative;
	unsigned int lo=1;
	unsigned int hi=n-1;
	while (hi>lo)
	{
		addFace(mesh, mesh.halfedge(hole[lo]).vertex_handle_, mesh.halfedge(hole[hi]).vertex_handle_, mesh.halfedge(hole[(hi+1) % n]).vertex_handle_);
		lo++;
	
		if (lo==hi) break;
		addFace(mesh, mesh.halfedge(hole[lo]).vertex_handle_, mesh.halfedge(hole[hi]).vertex_handle_, mesh.halfedge(hole[lo-1 ]).vertex_handle_);
		hi--;
	}

	return true;

//alternative:

// alternatively, use advancing front method

	double* tightness = new double [n];
	double* theta = new double [n];
	cout<<"Angles"<<endl;
	TriMesh::Point P;
	for (i=0; i<n; ++i)
	{
		tightness[i] = localDensity(mesh, mesh.halfedge(hole[i]).vertex_handle_);
		//theta[i] = angleBetweenPoints(mesh, mesh.halfedge(hole[ (n+i-1) % n]).vertex_handle_, mesh.halfedge(hole[i]).vertex_handle_, mesh.halfedge(hole[ (i+1) % n]).vertex_handle_);
		P = mesh.point(mesh.halfedge(hole[i]).vertex_handle_);
		theta[i] = mesh.calc_sector_angle(hole[i]);
		cout<<tightness[i]<<"  "<<theta[i]<<"  "<<P<<endl;
	}
	delete[] tightness;
	delete[] theta;
	return false;
} // bool closeHole(TriMesh& mesh, hehV& hole)

// http://stackoverflow.com/questions/5963954/fast-templated-c-octree-implementation

void testPolyTriangulation()
{
	LightVector<XY> points; // real asshole
	polygon poly;

/*	LightVector<XY> square;
	square.push_back(XY(-10, -10));
	square.push_back(XY(-10, +10));
	square.push_back(XY(+10, +10));
	square.push_back(XY(+10, -10));
//	square.push_back(XY(-20, 0));

	poly = square;

	poly.computeEdgesAngles();

	cout<<poly.inside(XY(0, 3))<<endl;
	cout<<poly.inside(XY(20, -5))<<endl;
*/
	
	points.push_back(XY(-42.02, 94.08));
	points.push_back(XY(-68.8, 76.87));
	points.push_back(XY(-97.31, 99.57));
	points.push_back(XY(-128, 106.17));
	points.push_back(XY(-151.17, 97.52));
	points.push_back(XY(-153.56, 76.98));
	points.push_back(XY(-148.75, 60.04));
	points.push_back(XY(-137.95, 43.05));
	points.push_back(XY(-107.17, 42.93));
	points.push_back(XY(-64.85, 41.39));
	points.push_back(XY(-27.59, 08.88));
	points.push_back(XY(-14.34, 43.03));
	points.push_back(XY(-2.66, 62.39));
	points.push_back(XY(-19.43, 80.02));

	poly = points;

	poly.computeEdgesAngles();
	poly.closeHole();

} // void testPolyTriangulation()


void printOff(triangulateio& tio)
{
	//std::ostream& io = cout;
	ofstream io;
	io.open("out.off");
	int i,j;

	io<<"OFF"<<std::endl;
	io<< tio.numberofpoints << "  " << tio.numberoftriangles << "   0" << std::endl;
	for (i=0; i<tio.numberofpoints; ++i)
	{
		io<<tio.pointlist[i * 2] << "   " <<tio.pointlist[i * 2+1] << "   0"<<std::endl;
	}

	io<<std::endl;
    for (i = 0; i < tio.numberoftriangles; i++)
	{
		io<<"3  ";
		for (j = 0; j < tio.numberofcorners; j++)
		{
          io<<tio.trianglelist[i * tio.numberofcorners + j]<< " ";
        }
		io<<std::endl;
	}
	io.close();
}

struct attributeType
{
	double dim3;
	int vhandle;
	OpenMesh::Vec3f normal;
};



void TriangleCall(TriMesh& mesh, hehV& hole)
{
// http://www.cs.cmu.edu/~quake/tripaper/triangle3.html

	struct triangulateio in, mid, out, vorout;
	int i,j;
	OpenMesh::VertexHandle vhandle;
	OpenMesh::FaceHandle fhandle;

	in.numberofpoints = hole.size();
	in.numberofpointattributes = (sizeof(attributeType) + sizeof(REAL) -1 )/ sizeof(REAL);

	in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL)); // points in 2 dimensions
	in.pointattributelist = (REAL *) malloc(in.numberofpoints * in.numberofpointattributes * sizeof(REAL));

	OpenMesh::Vec3f P;
	double sx(0), sx2(0), sy(0), sy2(0), sz(0), sz2(0);
	double stdx(0), stdy(0), stdz(0);

	for (i=0; i<in.numberofpoints; ++i)
	{
		vhandle = mesh.halfedge(hole[i]).vertex_handle_;
		P = mesh.point(vhandle );
		sx += P[0];
		sy += P[1];
		sz += P[2];

		sx2 += P[0]*P[0];
		sy2 += P[1]*P[1];
		sz2 += P[2]*P[2];
	}
	stdx += sx2/in.numberofpoints - sq(sx/in.numberofpoints);
	stdy += sy2/in.numberofpoints - sq(sy/in.numberofpoints);
	stdz += sz2/in.numberofpoints - sq(sz/in.numberofpoints);

	int dim1, dim2;
	if (stdx>stdy)
	{
		if (stdy>stdz)
		{ 
			dim1 = 0; dim2 = 1;
		} else 
		{
			dim1 = 0; dim2 = 2;	
		}
	} else
	{
		if (stdx>stdz)
		{ 
			dim1 = 0; dim2 = 1;
		} else 
		{
			dim1 = 1; dim2 = 2;	
		}
	}
	int dim3 = 3- dim1 - dim2;

	for (i=0; i<in.numberofpoints; ++i)
	{
		vhandle = mesh.halfedge(hole[i]).vertex_handle_;
		for (TriMesh::ConstVertexFaceIter vf_it=mesh.cvf_iter(vhandle); vf_it; ++vf_it)
			mesh.set_normal(vf_it.handle(), mesh.calc_face_normal(vf_it.handle()) );
	}

	// hole is projected on dim1 x dim2, this is where the stdev of the points is biggest
	TriMesh::Normal n;
	attributeType* iatr; 

	for (i=0; i<in.numberofpoints; ++i)
	{
		vhandle = mesh.halfedge(hole[i]).vertex_handle_;
		
		n = mesh.calc_vertex_normal(vhandle);
		cout<<n<<endl;

		P = mesh.point(vhandle);

		iatr = (attributeType*)(&in.pointattributelist[i*in.numberofpointattributes]);
		iatr->dim3 = dim3;
		iatr->vhandle = vhandle.idx();
		iatr->normal = n;

		in.pointlist[i*2] = P[dim1];
		in.pointlist[i*2+1] = P[dim2];
	}

	in.pointmarkerlist = NULL;

	in.segmentlist = (int*) malloc(in.numberofpoints * 2 * sizeof(int));
	in.segmentmarkerlist = (int*) malloc(in.numberofpoints * sizeof(int));
	for (i=0; i<in.numberofpoints; ++i)
	{
		in.segmentlist[i*2] = i;
		in.segmentlist[i*2+1] = i + 1;
		in.segmentmarkerlist[i] = 1;
	}
	in.segmentlist[in.numberofpoints*2-1] = 0;
	in.numberofsegments = in.numberofpoints;
 
	in.numberofholes = 0;
	in.numberofregions = 0;
	in.regionlist = NULL;
	in.numberoftriangles = 0;

//Make necessary initializations so that Triangle can return a triangulation in `mid' and a voronoi diagram in `vorout'.

	mid.pointlist = (REAL *) NULL;            // Not needed if -N switch used.
	// Not needed if -N switch used or number of point attributes is zero:
	mid.pointattributelist = (REAL *) NULL;
	mid.pointmarkerlist = (int *) NULL; // Not needed if -N or -B switch used.
	mid.trianglelist = (int *) NULL;          // Not needed if -E switch used.
	// Not needed if -E switch used or number of triangle attributes is zero: 
	mid.triangleattributelist = (REAL *) NULL;
	mid.neighborlist = (int *) NULL;         // Needed only if -n switch used.
	// Needed only if segments are output (-p or -c) and -P not used:
	mid.segmentlist = (int *) NULL;
	// Needed only if segments are output (-p or -c) and -P and -B not used:
	mid.segmentmarkerlist = (int *) NULL;
	mid.edgelist = (int *) NULL;         // Needed only if -e switch used.
	mid.edgemarkerlist = (int *) NULL;   // Needed if -e used and -B not used.

	vorout.pointlist = (REAL *) NULL;        /* Needed only if -v switch used. */
	/* Needed only if -v switch used and number of attributes is not zero: */
	vorout.pointattributelist = (REAL *) NULL;
	vorout.edgelist = (int *) NULL;          /* Needed only if -v switch used. */
	vorout.normlist = (REAL *) NULL;         /* Needed only if -v switch used. */
	
	triangulate("pzAevn", &in, &mid, &vorout);

	printf("Initial triangulation:\n\n");
	//  report(&mid, 0, 1, 1, 1, 1, 0);
	printOff(mid);

	// Attach area constraints to the triangles in preparation for refining the triangulation.

	// Needed only if -r and -a switches used:
	mid.trianglearealist = (REAL *) malloc(mid.numberoftriangles * sizeof(REAL));
	mid.trianglearealist[0] = 3.0;
	mid.trianglearealist[1] = 1.0;

	// Make necessary initializations so that Triangle can return a triangulation in `out'.

	out.pointlist = (REAL *) NULL;            // Not needed if -N switch used.
	// Not needed if -N switch used or number of attributes is zero:
	out.pointattributelist = (REAL *) NULL;
	out.trianglelist = (int *) NULL;          // Not needed if -E switch used.
	// Not needed if -E switch used or number of triangle attributes is zero:
	out.triangleattributelist = (REAL *) NULL;

	// Refine the triangulation according to the attached triangle area constraints.

	triangulate("prq22zBP", &mid, &out, (struct triangulateio *) NULL);

	printf("Refined triangulation:\n\n");
	//  report(&out, 0, 1, 0, 0, 0, 0);
	printOff(out);

	// add new points
	double w, dist, cord3;
	attributeType* matr = (attributeType*)(mid.pointattributelist);
	attributeType* oatr = (attributeType*)(out.pointattributelist);
	int nbNewPoints = out.numberofpoints - in.numberofpoints;
	int nbNewTriangles = out.numberoftriangles - in.numberoftriangles;
	LightVector<TriMesh::FaceHandle> newFaces(nbNewTriangles);
	LightVector<TriMesh::VertexHandle> newPoints(nbNewPoints);

	for (i = mid.numberofpoints; i<out.numberofpoints; ++i)
	{
		dist = 0;
		cord3 = 0;
		for (j=0; j<mid.numberofpoints; ++j)
		{
			w = 1.0 / (sq(mid.pointlist[2*j]-out.pointlist[2*i]) + sq(mid.pointlist[2*j+1]-out.pointlist[2*i+1]));
			dist += w;
			cord3 += w* matr[j].dim3;
		}
		cord3 /= dist;
		oatr[i].dim3 = cord3;
		//out.pointattributelist[2*i] = cord3;
		P[dim1] = out.pointlist[2*i];
		P[dim2] = out.pointlist[2*i+1];
		P[dim3] = cord3;
		vhandle = mesh.add_vertex(P);
		newPoints.push_back(vhandle);
		//out.pointattributelist[2*i+1] = vhandle.idx();
		oatr[i].vhandle = vhandle.idx();
	}
	// add new triangles
	int pointA, pointB, pointC;
	for (i = in.numberoftriangles; i < out.numberoftriangles; i++)
	{
		cout<<"Triangle"<< i << "  ";
		pointA = out.trianglelist[i * out.numberofcorners];
		pointB = out.trianglelist[i * out.numberofcorners + 1];
		pointC = out.trianglelist[i * out.numberofcorners + 2];
		cout<<pointA<<"  "<<pointB<<"  "<<pointC<<endl;

		fhandle = addFaceSimple(mesh, OpenMesh::VertexHandle(out.pointattributelist[2*pointA+1]),
							OpenMesh::VertexHandle(out.pointattributelist[2*pointB+1]),
							OpenMesh::VertexHandle(out.pointattributelist[2*pointC+1]));
		newFaces.push_back(fhandle);
	}
	cout<<"nbNewTriangles="<<nbNewTriangles<<endl;
	cout<<"nbNewPoints="<<nbNewPoints<<endl;
	// Free all allocated arrays, including those allocated by Triangle.
	free(in.pointlist);
	free(in.pointattributelist);
	free(in.pointmarkerlist);
	free(in.regionlist);
	free(mid.pointlist);
	free(mid.pointattributelist);
	free(mid.pointmarkerlist);
	free(mid.trianglelist);
	free(mid.triangleattributelist);
	free(mid.trianglearealist);
	free(mid.neighborlist);
	free(mid.segmentlist);
	free(mid.segmentmarkerlist);
	free(mid.edgelist);
	free(mid.edgemarkerlist);
	free(vorout.pointlist);
	free(vorout.pointattributelist);
	free(vorout.edgelist);
	free(vorout.normlist);
	free(out.pointlist);
	free(out.pointattributelist);
	free(out.trianglelist);
	free(out.triangleattributelist);
}





void main()
{
//	testPolyTriangulation();

	//tutorial1();
	//tutorial2("tetrahedron.off", "tetrahedron2.off",5);
	//tutorial1Hu();

//	Mat3 m(1,2,3,4,4,6,7,8,9);
//	Mat3 m2 = m.inverse();
//	cout<<m.determinant()<<endl;
//	cout<<m2;

	// TriMesh mesh = createSphere(2.0f,4); OpenMesh::IO::write_mesh(mesh, "tetraSphere.off");
	//mesh = createTetra(2); OpenMesh::IO::write_mesh(mesh, "tetra.off"); readmesh(mesh, "tetra.off");

	TriMesh mesh;
	readmesh(mesh, "lyukas.off");

	LightVector<hehV> holes;
	findHoles(mesh, holes);
	displayHoles(mesh, holes);

	// suboptimal, add normals to all faces while we need it only for the 1-ring around the hole
	if ( ! mesh.has_face_normals())
		mesh.request_face_normals();
//	mesh.update_normals();

	TriangleCall(mesh, holes[0]);
	TriangleCall(mesh, holes[1]);
	TriangleCall(mesh, holes[2]);
//	printOff(mesh, holes[2]);

	mesh.release_face_normals(); // always release after request

//	closeHole(mesh, holes[1]);
//	closeHole(mesh, holes[0]);
//	closeHole(mesh, holes[2]);
	holes.clear(); 
	
	findHoles(mesh, holes);
	displayHoles(mesh, holes);

/*	mesh.request_face_normals();
	if (mesh.has_vertex_normals() ) mesh.update_vertex_normals();

  for (; v_it!=v_end; ++v_it)
    this->set_normal(v_it.handle(), calc_vertex_normal(v_it.handle()));
*/



/*	readmesh(mesh, "icosahedron.off");
	OpenMesh::Vec3f P(-2,0.2,2);
	OpenMesh::FaceHandle fh = faceClosestToPoint(mesh,P);
	cout<<fh.idx();
	
	extendMesh(mesh, fh, P, true);
	OpenMesh::IO::write_mesh(mesh, "ico--.off");
	*/

	//cout<<"Volume = "<<meshVolume(mesh)<<endl;
	OpenMesh::IO::write_mesh(mesh, "lyukasvolt.off");
} // void main()

